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ABSTRAC T 


While Krylov and Bogolyubov used harmonic functions in their averaging 
method for the approximate solution of weak nonlinear differential equa- 
tions with oscillatory solution, we apply a similar averaging technique 
using Jacobi elliptic functions. These functions are also periodic and 
are exact solutions of strong nonlinear differential equations. The 
method is used to solve nonlinear differential equations with linear and 
nonlinear small dissipative terms and/ or with time dependent parame- 
ters. It is also shown that quite general dissipative terms can be 
transformed into time- dependent parameters. 

As a special example, the Langevin (collisional) equation of motion of 
electrons in a neutralizing ion background under the influence of a time 
and space-dependent electric field is presented. The method may also 
be used for nonlinear control theory, dynamic and parametric stabiliza- 
tion of nonlinear oscillations in plasma physics, etc. 


*On leave of absence Institute for Theoretical Physics, University of Innsbruck, 
Innsbruck, Austria. 


AVERAGING METHOD FOR THE SOLUTION 
OF NONLINEAR DIFFERENTIAL EQUATIONS 
WITH PERIODIC NONHARMONIC SOLUTIONS 


1. INTRODUCTION 

In these days, nonlinear differential equations are solved mainly numerically by 
using digital computers. This is an efficient and well-paved way to obtain a 
particular solution. However, if the system of differential equations is the 
Lagrange subsidiary system of characteristic equations for a first-order partial 
differential equation such as the Vlasov equation or another more sophisticated 
kinetic equation (in which the author was originally interested in connection with 
collisional Landau damping and the saturation of the two-stream instability) or 
if the influence of the parameters of the equation and/ or of the initial or boundary 
conditions on the solution are of interest, the computer solution can be very ex- 
pensive or even impossible if the computing time exceeds reasonable time 
intervals. 

Actually, some authors have recently proposed analytic methods for the solution 
of nonlinear, but nearly linear, differential equations. Krylov and Bogolyubov [ 1] , 
[ 2] investigated equations of the type 

x + co 2 x = e F(x, x) (1) 


where e is a small parameter. The methods starts from the so-called gener- 
ating solution 


x - A sin (cot + cp), x - A co cos (co t + cp) 


( 2 ) 


which satisfies (1) exactly in zero order (e = 0). 

In order to solve (1), it is assumed that the constants of integration A and <p de- 
pend on time, so that in (2) A - A (t), <p — cp (t). Expressing 1 (x,x) into a Fourier 
series in cp and assuming that the parameter e is small, so that amplitude A and 
phase qp change very slowly during one period of the oscillation; i.e., that 


A/A « co, q>/cp « co 


(3) 


one obtains in first order of e by averaging over one period 


1 


(4) 



e 

O) 


£ 

2 rr 



sin qp, A a) cos qp) cos qp dcp 



_e_ J_ 
Ago 2tt 



sin qp, A a; cos (?) sin qp dcp 


(5) 


where A and <p are assumed to be time independent under the integrals. Also 
higher order solutions can be obtained [ 22] . 

This method which is, however, restricted to equations of the type (1) (i.e., to 
nearly linear equations) has been used extensively in plasma physics, theory of 
oscillations, control theory, etc., [3] to [ 7] . The method was also used to solve 
partial differential equations, [3] and [4] . 

Kruskal ( 8] extended the Krylov -Bogolyubov method to solve equations of the 
type 


or 


F(x, € ) 


no 



n = 0 


d n F (x, c) 
de n 


e = 0 


( 6 ) 


x - F(x, x, e) 


The solutions of these fully nonlinear equations are based on recurrent solutions 
and are given in the form of power series of the smallness parameter e . 

In this paper , we are investigating an averaging method for the solution of equa- 
tions of the type 


x + co 2 f (x) - e F(x, x) (7) 

using elliptic functions. We are starting with the exact solution of the fully non- 
linear equation (unperturbed equation) 

x + CO 2 f (x) = 0 (8) 

without developing f (x) into a power series with respect to the smallness 


2 


parameter e. Actually, one could bring (7) into the form (8) and apply Kruskal's 
method. This paper may therefore be considered as a special example of 
Kruskal’s theory, using not a transformation technique but a direct Krylov- 
Bogolyubov type technique, averaging over one period of elliptic functions and 
giving concrete examples. Furthermore, a generalization to time dependent 
parameters is given. 


2. ELLIPTIC FUNCTIONS AS PERIODIC NONHARMONIC SOLUTIONS 
OF NONLINEAR DIFFERENTIAL EQUATIONS 

The Krylov-Bogolyubov technique is based on harmonic solutions. If, e.g., one 
considers Duffing's equation 


X + OJ 2 X = € X 3 


( 9 ) 


or Einstein's equation for the perihelion shift 

x + a ; 2 x = e x 2 + a (10) 


then the Krylov-Bogolyubov technique starts from (2) and is valid only for e < < 1. 
One obtains from (4) the result A = const and the frequency modification ("ampli- 
tude dispersion") is contained in (5). It is, however, possible to solve (9) exactly 
and any other equation of type (8) for any e . By multiplying (8) by x and inte- 
grating twice, one obtains 


t - t, 


■ f 


dx 




0 f (x) dx - 2E 


( 11 ) 


where t 0 and E are integration constants. If f (x) is a polynomial of degree up to 
three or a simple harmonic function like sin £x, then (11) is an elliptic Integral 
and its inverse function may be expressed by a Jacobi elliptic function [9] to [12]. 

We are now going to consider the Langevin equation of motion of electrons in a 
periodic space-dependent electric field [ 13] . This equation is of importance in 
plasma physics, in kinetic theory (as Lagrange characteristic equation of colll- 
sional kinetic equations), etc., but also as a form of the Froude equation for 
rolling ships or the damped pendulum. It reads (o>, 1, e are given constants) 


3 



x + co 2 • sin £x 


e F (x, x) 


( 12 ) 


where £ is a collision frequency (e.g., electrons in a neutralizing ion background). 
The solution for e = 0 (unperturbed equation, generating solution in the sense of 
Krylov- Bogolyubov) is 


2 . 

x(t) ~ arcsin [k sn(yEv//,k)] (13) 

where sn = the Jacobi elliptic sine function, \p - cot + <p and k, the modulus of 
the elliptic function (amplitude of the oscillation), and q> are integration constants. 


3. AVERAGING METHOD WITH ELLIPTIC FUNCTIONS 

In order to solve (12), we now replace according to Krylov-Bogolyubov [ 1] 

k - k (t), qp - cp(t) (14) 


and set up as generating solution (13) and 


x 


2k a 

rr 


cn(Vf qy ,k) 


(15) 


which may be obtained by deriving (13) with k and cp kept constant and then ap- 
plying (14). Also the relations 


1 - k 2 sn 2 - di 2 , = I? cn dn (16) 

were used [9] . Differentiating now (13) and observing (14), then equating it to 
(15) gives 

, dsn 

k m * k 'ST 

" ' E yr cn dn < 17 > 


Since <p = qp (k) for Jacobi functions [ 9] , we have to use the total derivative with 
respect to k when using d/ dt. Equation (17) is independent of the form of F (x, x). 
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Differentiating (15) and substituting into (12), we obtain 


kcn-k VFsndn* <p + kk 


d cn 
dk 


eVf -x 

“ ~2Z F(X,X) 


(18) 


Here we used 


B cn 

B \p 


VF sndn and sin £ x 


£ x £ x 
2 sin ~2 cos - 


2ksncos arcsin [] 


2ksncos arccos 


V'i - n 2 


(19) 


and (16). 

Solving (17) and (18) for k and , we receive 

£ = F ( x -’ c ) cn(VF^k) (20) 


Here use has been made of 


sn 2 + cn 2 


1, 


sn 


d sn 
~dk~ 


d sn 2 
dk 2 


Furthermore, we have 


<p - 


€ F ( x, x) 
2 co dn 


sn + k 


d sn 
dk 


( 21 ) 


( 22 ) 


Since the Jacobi functions sn, cn, and dn are periodic with the period 4K where 
the quarter period [ 9] 


K(k) = 


r — 

l 1 FT- 


dd 


Y 1 - k 2 sin 2 0 


(23) 
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is the complete elliptic integral of the first kind, we may average (20) and (21) 
without any Fourier series expansion. Defining u =V7 cp and 

1 f 4K 

<"> ~ 4K J Q ’ * d ^ ’ ( 23 ) 


we now calculate <ic> and <q?> in analogous manner to (4) and (5). We then have 
from (20) 

r4K 

F(snu,cnu) cnudu (24) 

o 


where again <p and k are considered to be constant under the integral. Further- 
more, from (22) 




r4K 


£ J_ 

oj 8K 


sn u + k 


F (sn u, cnu)- 


dsn 


du 


in u 




(25) 


We consider now some special cases, 
a. F (x, x) = F (x) = F (sn u) 

This case is of no interest, since (12) is then of the form (8) and can be integrated 
exactly. (See (11).) From (24), we actually get k = 0, since cn u du = d sn u/dn u 
or with sn u = y, using (16) we have 



4 > (m u) 


4K 

0 


4 > (m 4K) - <£(m0) 


0 (26) 


b. F (x, x) = F (x) = F (cn u) 

We then have from (25) using cn u = y and 


dn 2 


( 1 - k 2 ) k 2 cn 2 . 


du 
dn u 


d cn u 
sn u dn 2 u 


(27) 
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so that the theorem by Krylov and Bogolyubov that dissipative terms in first 
order do not modify the phase (frequency) is also valid here. 


4. SOME SPECIAL EXAMPLES 
We now present some applications of the method, 
a. Linear Damping 
We have 


F (x) = -x = - 


2ku 


cn(u, k) 


(29) 


Then <$> = 0 from (28) and (25) , (24) give (9] 


/l d In 

Vk K < k > St 


4 yr Jo 


r 4K e 1 

J cn 2 u du = - yr * — [E(k) - ( 1 - k 2 /K(k)] (30) 


where 


r w/2 

E(k) = j yi - k 2 sn 2 d d 


v 


(31) 


is the complete elliptic integral of the second kind. Using the identity (9) 


kK = rk [E(k) - (1 - k 2 )K(k)] * H’(k) 


(32) 


we then have 


vr 


(t - t 0 > = 


j¥(kf dk = ln £ E < k ) - (1 - k2 ) K(k)3 (33) 


which gives k = k (t). 


b. Quadratic Damping 


We have 


F(x) = -x 2 


4 k 2 u> 2 
£ cn 


(34) 


and from (24) we get k = 0. This is understandable because, from physical con- 
siderations, the damping function should be an odd function of x; e.g., F(x) = ixlx. 
See also [ 2] or F = x + x 2 . However, variable damping is of more interest. 

c. Van der Pol Damping 

We generalize the Van der Pol equation 

x + x- e(l - x 2 ) x = 0 (35) 

to the form (12); l.e., 

x + co 2 sin ? x = — 6 ( x 2 — 1 ) x * f F (36) 

This would represent, for example, a bounded inhomogeneous plasma model. See 
Lashinsky [3] . 

From (24), we again get exactly (30) and (33). 

A warning might be useful: Before applying the method described here, one has 
to determine if (8) has a periodic solution at all. There are cases with "big" e 
in which (8) has a periodic solution, but (7) does not have. An example is equa- 
tion [ 14] 


x-ax+cx 3 = -ex (36) 

For e = 0, an elliptic function is the solution. For e 2 < 8 a, we have a stable 
focus and a damped oscillation and the method of Section 3 can be applied. For 
e 2 > 8 a, we have a stable node and no oscillation at all. 
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5. GENERALIZATION FOR NONE LLIPTIC PERIODIC FUNCTIONS 
AND FOR TIME-DEPENDENT PARAMETERS 

If f (x) in (7) is such that the solution of (8) is not an elliptic function (it could be 
a hyperelliptic or an Abelian function with more than four periods) , then a similar 
method can be devised if, and only if, x(t) is periodic. 

In some other cases, the parameter co in (12) is time dependent; e.g., in the in- 
vestigation of nonlinear Landau damping [ 15} and [ 16] . On the other hand, a 
dissipative equation of the type (7) in the special form x 

x + co 2 f (x) = e g(x) h(x) (37) 

(here & is constant) with generalized Rayleigh damping 

g(x) - a x 3 + f3 x 2 + y x + 8 (38) 


and arbitrary h (x), x = x ( f ) , £ (V) can be transformed (see Appendix) into 


£ + — ■ P(£) = 0 (39) 

nr * 


so that even nonlinear and variable dissipation terms and therefore in special 
cases (7) and (12) can also be reduced to a "dissipationless" nonlinear oscilla- 
tion equation with time-dependent parameter of the form 

x + to 2 ( t ) H (x) = 0. (40) 


It does not appear to be possible to solve this equation exactly by quadratures 
[23]. 

f (x) of (37) is mow no more restricted to a polynomial of 3rd degree and co (t) is 
given either by the physical problem (e.g., [15] to [ 16] ) or from the transforma- 
tion (e.g., co = t _1 ). 

If, however, co - o;(t) is given and dissipation is present (e.g., in the Langevin 
equation, see [13] , [17] and Section 6), then we have instead of (7), (12) resp 
(40), the more general equation 
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X + (j?- (t) f (x) 


e F (x, x) 


(41) 


which cannot be transformed into (40) . (See Appendix.) 

hi order to solve (41), we first solve (8) for oo = const. Let the solution be 

(42) 

(43) 


where 


x(t) - y(0,k) 


0(0 


Cl) t + CP 


and where <p and k are integration constants. Using (14), we then set up as gen- 
erating solution 


x(t) - y (0 t, k(t) , k(t)) (44) 

x(t) - o;(t) = c<j(t) • y^ (45) 


Since the argument 0 of Jacobi functions depends on k, we have to assume 
3 cp / 3 k ^ 0 also in the more general case. Equating (45) with dx / dt of (44) 
gives 0 analogous to (17) . Differentiation of (45) and substitution into (41) gives 
k. Solving for 0 and k gives 


. Sk (EF ~ y *“> . 

* = ijN 


(46) 


>>(£ F - y f'u) 
oN 


(47) 


where 


N - yk ~ y* y*k 


(48) 


Using the identity 
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dN 

d i// 


(49) 


d y d 

d k ^'P d k 

and differentiating y ^ + f = 0 with respect to 0 , respective with respect to k 
and subtracting, one may show that dN/d0 = 0 (i.e., N is independent of 0, so 
that it can be put before the integral J . . d 0 .) 

We now assume that the time-dependent parameter co varies slowly with time in 
(41). For (40), this assumption cannot be made if a ; (t) stems from dissipative 
terms. As a measure of this slowness, we introduce another small parameter n 
and define [18] as a stretched time variable 


# = M, 


d CD 

d t 


d Cl) 

^ dd 


(50) 


so that two scales defined by e and , u are now involved. 
Assuming 


0(0) % 0(T), l « ~ 



1 d Cl! 1 

« d# <<: 7 Tt 


(51) 


where T is the smallest period of the periodic function (42), we now average (46) 
and (47) over one period assuming that <p, 0, k, co remain constant during T. We 
then have 


< q>> 


e f T dy 
wNT J, dk 


F (y, y) d0, 


< 0 > 


co + 


(52) 


<k> 


€ 

co NT 



Cl) 

CD NT 



(53) 


which are the equivalents of (25) and (24). In the derivation, the relations 


d y - y^ d 0 , t 


0 - <P 

CD 
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where cp , are constant under the integral and 


f T d_y 

Jo dk 


dy 


d r T 

^ V 


y d y 


_d_ 
d k 


[H(y)lJ 


0 


were used. For F(y,y) - F (y), the first r.h.s. term in (53) vanishes. For 
F (y, y) = F (y), no general conclusion could be reached. 

Sometimes the integrals in (52) and (53) are difficult to obtain. It might then be 
useful to start with the differential equations for the function t(x), which is the 
inverse function of the function x(t) solving (41). Before the averaging process, 
it is necessary to return to the original functions. 


6. THE LANGEVIN EQUATION OF MOTION 

We are now applying the method described on the Langevin equation of motion of 
an electron. The equation reads [ 13] 

x + <^i 2 (t) sin f x = -ex (54) 

Here, oj (t) is given from either the Maxwell equations or from the energy trans- 
fer rate between electrons and the electromagnetic field, [ 15] to ( 16] . Using 
(13) and (29) as generating solution, we get from (52) and (53) the equations 


<<p> - CO 

(55) 

|e + — j ken 2 (]/i \p , k) 

(56) 


so that 


w(t) [E(k) - K(k) (1 - k 2 )] - const. e“ e * (58) 

determines now k(t). With this result, Denavit's calculations on the collisionless 
Landau damping [6] could be extended to collisional damping [ 19] to [20] . 
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APPENDIX 


We are presenting here the transformation of (37) into (39). The method is due 
to Abel and may be found in Kamke [ 21] . The consecutive transformations are 

x(t) = u(t) = u(x), X = u’(x)*u(x), u’ = u = V (59) 

which give 

-v’ + v 3 a; 2 f (x) = eg[~Jh(x)v 3 (60) 


where g(v) is given by (38), Equation (60) is an Abelian equation. Using 
V (x) = W(x) • Tj(f), £(x) = - / w(x) e y h(x) dx , 

(61) 

w(x) = exp [- J e /3h(x) dx] , 


we obtain for a - 0 (for a ^ 0 see Kamke [ 21] for a more general transforma- 
tion into another equation than (39) ) the Abel equation in standard form 


dr} 

d? 


-^P(£) + r l 2 


(62) 


where P (x(£) ) is given by 


a) 1 f (x) -eh h(x) r 

P(x) = — ■ ■ T - 7 "t — — e ' £ ^ h(x > d * 

v ' e y h(x) 


(63) 


x(£) is given by (61). Using then 


d£ 

df E £’< r > 


1 


T l(€) 


(64) 


which defines r, we obtain (39). 
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